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Abstract 

Constitutive equations are derived for the time-dependent behavior of particle- 
reinforced elastomers at isothermal loading with finite strains. A rubbery polymer 
is modelled as a network of macromolecules bridged by junctions, where the junc- 
tions can slip with respect to the bulk material under straining. A filled rubbery 
compound is thought of as an ensemble of meso-regions where sliding occurs with 
different rates. Stress-strain relations for a particle-reinforced rubber are developed 
by using the laws of thermodynamics. For uniaxial tensile relaxation tests, these 
equations are determined by three adjustable parameters. To find the experimental 
constants, three series of tests are performed at longitudinal strains in the range 
from 100 to 250 %. By fitting observations, the effects of pre-loading and ther- 
mal recovery are analyzed on the nonlinear viscoelastic response of natural rubber 
reinforced with carbon black. 
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1 Introduction 



This paper is concerned with modeUing the time-dependent behavior of particle-reinforced 
elastomers at isothermal loading with finite strains. The viscoelastic and viscoplastic re- 
sponses of filled rubbers have been a focus of attention in the past decade, which may 
be explained by numerous applications of rubber-like materials in industry (vehicle tires, 
shock absorbers, earthquake bearings, seals, flexible joints, solid propellants, etc.), see, 
e.g., Govendjee and Simo (1992), Johnson and Stacer (1993), Hausler and Sayir (1995), 
Johnson ct al. (1995), Aksel and Hiibner (1996), Holzapfel and Simo (1996), Lion (1996, 
1997, 1998), Spathis (1997), Bcrgstrom and Boycc (1998), Ha and Schapcry (1998), 
Kaliskc and Rothert (1998), Reese and Govindjee (1998), Septanika and Ernst (1998), 
Miehe and Keck (2000), Wu and Liechti (2000), Haupt and Sedlan (2001), to mention a 
few. Although the number of publications on this subject is rather large, some important 
issues remain, however, obscure. 

One of the main questions to be resolved deals with the physical nature of the time- 
dependent response of rubber-like materials. It is worth noting that there is no agreement 
in the literature about the micro-mechanisms for stress relaxation in unfilled elastomers. 

The viscoelastic behavior of rubbery polymers is conventionally modelled within the 
concept of temporary networks (Green and Tobolsky, 1946). A polymer is treated an a 
network of macromolecules bridged by junctions (chemical and physical cross-links, en- 
tanglements and inter-chain links induced by van der Waals forces). With reference to the 
affinity hypothesis, junctions are assumed to be rigidly connected with the bulk medium, 
whereas the time-dependent response is associated with breakage of active strands (whose 
ends are linked to separate junctions) and reformation of dangling strands (where stresses 
totally relax before these chains merge with the network). 

Another scenario for stress relaxation in rubbery polymers is based of the concept 
of non-affine junctions that presumes the network of macromolecules to be permanent 
and explains the time-dependent response at the macro-level by slippage of junctions 
with respect to the bulk material at the micro-level. This approach was first suggested 
by Phan-Thien and Tanner (1977) and Phan-Thien (1978) for polymer melts and was 
applied by Desmorat and Cantournet (2001) to study the behavior of filled elastomers. 

A version of the concept of non-affine networks, where junctions do not slide, but 
polymeric chains are allowed to slip with respect to entanglements was proposed by Miehe 
and co-authors for the analysis of stresses in rubbery vulcanizates (Michc, 1995; Michc 
and Keck, 2000, Lulei and Michc, 2001). A similar model (where chains slip with respect 
to filler particles instead of entanglements) was recently introduced by Hansen (2000) 
with reference to results of molecular dynamics simulations and atomic force microscopy 
measurements. 

The above theories concentrate on relatively slow relaxation processes when stresses 
relax at the time-scale of structural rearrangement in the entire macromolecule. A rel- 
atively fast relaxation process (whose characteristic time is associated with the time for 
rearrangement of a mechanical segment in a chain) was studied by Drozdov (2001) by 
using a micro-mechanical scenario for the time-dependent response of rubbery polymers 
similar to the tube model for concentrated polymeric solutions and melts (Doi and Ed- 
wards, 1986). The difference between the two approaches is that the conventional tube 
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theory is based on the reptation concept for polymeric chains (this kind of diffusive mo- 
tion is prevented by cross-hnks in rubbery vulcanizates) , whereas the concept of wavy 
tubes ascribes the viscoelastic behavior of rubbers at the macro-level to mechanically 
induced changes in the concentration of regions with high molecular mobility in an indi- 
vidual strand at the micro-level. A transition mechanism from relaxation processes at the 
length-scale of a statistical segment to those at the length-scale of a chain in a polymer 
melt was recently discussed by Herman (2001). 

This brief discussion of basic concepts in finite viscoelasticity and viscoplasticity of 
rubbcr-likc materials demonstrates that our knowledge of micro-mechanisms for the time- 
dependent response of unfilled elastomers is far away from being exhausted. The physics 
of the viscoelastic and viscoplastic behavior of particle-reinforced rubbers remains rather 
unclear. This conclusion is confirmed by a number of micro-mechanical theories that have 
been proposed in the past decade to describe the viscoelastic response of filled elastomers. 

On the one hand, Aksel and Hiibner (1996), Kaliske and Rothcrt (1998) and Hansen 
(2000) explained the time-dependent behavior of particle-reinforced rubbers by sliding of 
polymeric chains along and their dc- wetting from the surfaces of filler clusters. This con- 
cept goes back to the model proposed by Dannenberg (1975), which, in turn, is equivalent 
to the neglect of the ability of a host matrix to relax, the hypothesis suggested about half 
a century ago by Bucchc (1960, 1961) to explain the Mullins effect. 

On the other hand, Liang et al. (1999), Clarke et al. (2000) and Leblanc and Car- 
tault (2001) evidenced that the effect of particles and their aggregates on the viscoelastic 
response of rubbery compounds is not crucial. Leblanc and Cartault (2001) performed dy- 
namic shear tests on uncured styrene-butadiene rubber reinforced with various amounts 
of carbon black and silica and reported rather limited time-dependent responses of the 
compounds. These observations allowed them to infer that inter-molecular interactions 
play the key role in the viscoelastic behavior, whereas the influence of flller is of second 
importance. Liang et al. (1999) reported experimental data in dynamic tensile tests 
on ternary composite of polypropylene and EPDM (ethylene propylene diene monomers) 
elastomer filled with glass beads. They demonstrated that treatment of beads by silane 
coupling agent did not affect mechanical damping and that the loss modulus of the com- 
pound with untreated beads decreased with the concentration of filler (in contrast with 
the de- wetting concept). Clarke et al. (2000) carried out tensile relaxation tests on two 
elastomers (natural polyisoprene rubber and synthetic polyacrylate rubber) and revealed 
significant relaxation of stresses in the unfilled natural rubber. 

Our experimental data in tensile relaxation tests on natural rubber reinforced with 
high abrasion furnace black demonstrated comparable decreases in the longitudinal stress 
for a filled [45 phr (parts per hundred parts of rubber) of carbon black (CB)] and an 
unfilled [1 phr of CB] compounds (Drozdov and Dorfmann, 2001). Because the amount of 
stress relaxed during 1 hour at room temperature in the filled specimens exceeded that for 
the unfilled samples (approximately by twice), we conclude that a micro- mechanism for 
the viscoelastic response of filled elastomers is in between the above two extreme concepts: 
the relaxation process in particle-reinforced rubbers may be associated with rearrangement 
of macromolecules in the matrix, but the rearrangement process is strongly affected by 
local inhomogeneities induced by the presence of filler. 

In accord with this picture, we model a filled rubber as a network of macromolecules 
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bridged by junctions. To account for the effect of filler particles and their agglomerates 
on the time-dependent response of the network in a tractable way, the network is treated 
as an ensemble of meso-domains, where rearrangement of strands occurs with different 
(mutually independent) rates. The characteristic size of these regions (roughly estimated 
as several hundreds of nm) is assumed to substantially exceed the radius of gyration for 
macromolecules and the average radius of particles, on the one hand, and to be noticeably 
less than the dimensions of a macro-specimen. The ensemble of meso-domains reflects 
local inhomogeneities in the rubber-like material caused by impurities, non-homogeneity 
in the spatial distribution of a cross-linker and a filler, segregation of short chains to 
interfaces between the filler clusters and the host matrix (Carlier et al., 2001), local 
fluctuations of elastic moduli associated with the presence of hard cores around flUer 
particles and their clusters, as well as the enhancement of the rearrangement process 
for strands in the close vicinities of the flUer aggregates (Karasek and Sumita, 1996). 
The work is conflned to flUed rubbers with the volume fractions of particles, below 
the percolation threshold, which implies that phenomena associated with aggregation of 
clusters into an inflnite secondary network (Yatsuyanagi et al., 2001) are not taken into 
account. 

With reference to the concept of sliding junctions, different meso-regions are charac- 
terized by different potential energies, u, for slippage of junctions with respect to the bulk 
material. The rate of slippage, F, is expressed in terms of the potential energy, cj, by the 
Boltzmann formula. This implies that the sliding process is entirely determined by the 
attempt rate, Fq, and the distribution function, p{ijj), of potential energies for sliding. 

To account for the influence of strain intensity on the rate of stress relaxation, we 
distinguish between active meso-domains, where junctions can freely slide, and passive 
meso-regions, where sliding is prevented by surrounding macromolecules (because local 
concentrations of cross-links and van der Waals links between strands in these domains 
exceed their average values). Under stretching, some van der Waals links between strands 
break in passive meso-domains, which implies that these regions become active. The 
ability of a passive meso-region to be activated is determined by a measure of damage, 
a, that determines the relative number of extra links between macromolecules which do 
not permit the sliding process to start. The present study focuses on loading programs in 
which passive meso-regions are transformed into the active state when the stress intensity 
grows. Healing of inter-chain links that results in the inverse transformation at unloading 
is beyond the scope of this paper. 

The objective of this paper is three-fold: 

1. To develop constitutive equations for the time-dependent response of a particle- 
reinforced elastomer modelled as an ensemble of meso-regions with various potential 
energies for sliding of junctions. 

2. To report experimental data in three series of tensile relaxation tests on CB flUed 
natural rubber at various longitudinal strains in the range from 100 to 250 %. 
The first series of tests was carried out on virgin specimens, in the second series 
the samples were initially pre-loaded, whereas the third series dealt with the same 
specimens after thermal recovery at an elevated temperature. 
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3. To find adjustable parameters in the stress-strain relations by fitting observations 
and to rationalize the influence of pre-loading and thermal recovery on the nonlinear 
viscoelastic behavior of filled elastomers in terms of the constitutive model. 

The study is organized as follows. Sliding of junctions in active meso-regions and 
transition of passive meso-domains into the active state are discussed in Section 2. Kine- 
matic equations for slippage of junctions with respect to the bulk medium arc developed 
is Section 3. The mechanical energy of a particle-reinforced elastomer is determined in 
Section 4. Stress-strain relations are derived in Section 5 by using the laws of thermody- 
namics. Uniaxial extension of a specimen is analyzed in Section 6. Section 7 deals with 
the description of the experimental procedure. Adjustable parameters in the constitutive 
cqTiations are determined by fitting observations in Section 8, where a physical interpre- 
tation is provided for the findings. Some concluding remarks are formulated in Section 
9. 



2 Deformation of a filled elastomer 

A particle-reinforced elastomer is modelled as an ensemble of meso-domains (with ar- 
bitrary shapes) where junctions can slip with respect to the bulk material. Sliding of 
junctions is assumed to occur independently in different meso-regions. 

A mcso-domain is characterized by its potential energy cj, which is defined as follows. 
Slippage of junctions with respect to their positions in the bulk material is thought of as 
a viscous flow with the constitutive equation 

a = 7]e, (1) 

where a is the stress, e is the strain rate, and 77 > is some viscosity. The parameter 77 
in Eq. (1) is given by 

where G stands for an clastic modulus, and F is the rate of sliding. For definiteness, 
we suppose that G coincides for all meso-domains, whereas F accepts different values for 
different regions. The Boltzmann formula is adopted for the rate of sliding 

where Co is the potential energy of a meso-region, T is the absolute temperature, kB 
is Boltzmann's constant, and the pre-factor Fq weakly depends on temperature. The 
attempt rate Fq equals the rate of sliding for junctions at an elevated temperature {T — > 
00). Introducing the dimensionless potential energy. 



where Tq is some reference temperature, and disregarding the effect of small temperature 
increments, AT — T — Tq, on the sliding process, we find from Eq. (3) that 

F = Foexp(-u;). (4) 
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Equation (4) expresses the rate of slippage of junctions with respect to the bulk medium, 
r, in a meso- region in terms of its potential energy, lo. 

The potential energy lo characterizes the rate of sliding of junctions in a meso-domain 
where the sliding process is not affected by inter-chain interactions. The neglect of in- 
teraction between polymeric chains is a conventional hypothesis in the entropic theory 
of rubber elasticity (Treloar, 1975), where the number of available configurations for a 
chain is determined under the assumption that surrounding macromolecules do not re- 
strict the chain geometry. For a filled elastomer, interactions between polymeric chains 
may be disregarded only for a part of the ensemble of meso-domains. In the other part of 
meso-regions these interactions prevent sliding of junctions imtil micro-strains become suf- 
ficiently large to break inter-chain links (these links reflect van der Waals' forces between 
strands) that do not allow the sliding process to start. 

With reference to this picture, two groups of meso-domains in an ensemble are distin- 
guished: 

1. active regions where sliding of junctions is unrestricted, 

2. passive domains where slippage of junctions is inhibited by inter-chain interactions. 

The sliding process in an active meso-region is entirely determined by the potential energy, 
a;, of this domain. To characterize the ability of a passive meso-domain to become active, 
we introduce an additional variable, a, which is thought of as a measure for mechanically 
induced damage of inter-chain interactions. 

Any interaction between polymeric chains that prevents slippage of junctions with re- 
spect to the bulk material may be thought of as a temporary link between macromolecules. 
The weakening of interaction between strands driven by mechanical factors is modelled 
as rupture of temporary links under loading. 

Assuming the number of meso-regions with potential energy a; per unit mass of an 
elastomer to be rather large, we denote by (M(a;)) the average number of inter-chain 
links per unit domain before loading, and by M(i, to) the number of these links in a meso- 
domain in the deformed state at time t (the instant t = corresponds to the time when 
external forces are applied to a specimen). The measure of damage, a, is defined as the 
ratio of these quantities. 



Unlike conventional approaches to fracture of elastomers, where the initial measure of 

damage is assumed to vanish, whereas fracture of a macro-sample is associated with some 
positive value of the damage variable (conventionally, unity), Eq. (5) implies that a is 
non-negative in the stress-free state, and it decays to zero when a passive meso-region is 
transformed into an active one. 

It follows from Eq. (5) that the initial value of the damage variable, a, may be either 
higher or lower than 1 in passive meso-domains, whereas the average value of Q;(0,a;) 
equals unity. The condition q;(0,c<j) < 1 means that a passive meso-domain is initially 
damaged. This damage is associated with breakage of inter-chain interactions in meso- 
regions during the preparation process (at the stage of milling of a rubbery compound). 
On the other hand, the inequality q;(0, a;) > 1 implies that the number of links between 



a{i^ uS) 




(5) 
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macromolecules in a passive meso-domain in a virgin specimen exceeds its average value, 
which may be ascribed to local inhomogeneities in the distribution of a cross-linker. 

In accord with this picture, the distribution of passive meso-domains is described by 
the function Sp(t, a, uj) that equals the current number of passive meso-regions (per unit 
mass) with the potential energy uj where the damage variable equals a. Summing the 
number of meso-domains with various measures of damage, a, we find the concentration 
of passive regions with potential energy a;, 

/■oo 

Xp{t,uj)= Ep{t,a,uj)da. (6) 

The number of passive domains per unit mass of an elastomer, A^p, equals the sum of the 
numbers of passive meso-regions with various potential energies, uj, 

Np{t) = / Xp{t,u;)dw. (7) 

Jo 

Denote by Xa(t, u) the current number of active meso-domains with potential energy u. 
The quantities X^^ and Xp are connected by the conservation law 

X,{t,u;)+Xp{t,u;)^X{u;), (8) 

where X{uj) is the concentration of meso-regions with potential energy uj. Differentiation 
of Eq. (8) with respect to time implies that 

^{t,u;) = ^{t,u;)^-^{t,u;), (9) 

where '-f{t,uj) is the rate of transformation of passive meso-domains into active ones. 
Integrating Eq. (8) over cu and using Eq. (7), we find that 

N,{t) + Np{t)^N, (10) 

where 



iVa(t) = / x,{t,oj)dw (ii; 

Jo 

is the current concentration of active meso-domains and 

/■oo 

N= / X(uj)duj 
Jo 

is the total number of meso-regions per unit mass. 



3 Kinematic relations 

Let Tq be the radius vector of an arbitrary junction in the reference state, r{t, tq) its radius 
vector in the deformed state at time t, and r-a{t, Fq) the radius vector of this point in the 
unloaded state. In what follows, the argument Vq of the functions r and is omitted for 
the sake of simplicity. 
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The unloaded state of a mcso-region characterizes changes in the positions of junctions 
at time t driven of their shppage with respect to the bulk material (i.e., by slippage with 
respect to the initial configuration of a specimen) . The displacement vector for transition 
of a junction from its reference state to the unloaded state is given by 

Uu(t) = ru(t) - ro. 

Denote by Vor(t) the deformation gradient for transition from the reference state to the 
actual state, by Voru(^) the deformation gradient for transition from the reference state 
to the unloaded state, and by Vu(t)r(t) the deformation gradient for transition from the 
unloaded state to the deformed state. Here Vq and Vu(i) are the gradient operators in 
the initial and unloaded configurations, respectively. With reference to Drozdov (1996), 
the "left" gradient operators are employed: the image, dr{t), in the actual state of an 
infinitesimal vector, dvo, in the reference state reads 

dr{t) ^dro-Wor{t), 

where the dot stands for inner product. The deformation gradients are connected by the 
equation 

Vor(t) = Vor„(t)-Vu(t)r(t), (12) 

which expresses the chain rule for differentiation of the vector function r = r(t, Tq) with 
respect to Fq. 

The right Cauchy deformation tensor, C(t), for transition from the initial state to the 
deformed state is determined by the formula 



C(t) = Vor(t)- Vor(t) 



(13) 



where T denotes transpose. By analogy with Eq. (13), we introduce the right Cauchy 
deformation tensor for transition from the reference state to the unloaded state. 



Cu(t) = Voru(t) ■ Voru(t) 



(14) 



and the right Cauchy deformation tensor for transition from the unloaded state to the 
actual state. 



Ce{t) = Vu(t)r(i) • [Vu(i)r(i) 
It follows from Eqs. (12), (13) and (15) that 



(15) 



Ce(t) 



Voru(t) •Vor(t)- Vor(t) ■ Vor^{t) 



-T 



1 -1 



Vor,(i) •C(t)- Voru(i) 



(16) 



The left Cauchy deformation tensor for transition from the reference state to the deformed 
state reads 



B(i)= Vor(i) •Vor(i), 



(17) 
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whereas the left Cauchy deformation tensor for transition from the unloaded state to the 
actual state is given by 

Be(t) = [Vu(t)r(t)]^ ■ Vu(t)r(t). (18) 
Equations (12), (14) and (18) imply that 



Be(t) = [Vor(t)l^ ■ [Voru(t)l ^ ■ [Voru(t)l ' • Vor(t) 



Vor(t)l ■C-\t)-Vor{t). 



(19) 



Our purpose now is to determine the derivatives of the principal invariants, {k — 1,2), 
of the tensor Ce{t) with respect to time. Differentiation of Eq. (16) results in 



~~dt 



it) 



di 

+ 

+ 



Voru(t) -Cit)- Vor,(i) 



-T 



Bearing in mind that 



„ , ,1-1 dC, , r„ , ,1-T 

Vor,(t)J ■ —{t) ■ [Vor,{t)\ 
Voru(i)]"' • C{t) ■ or n{t)Y^ . 



^(t)=v.(t). 



(20) 



where Vu(t) is the velocity vector for sliding of junctions, and using the chain rule for 
differentiation, we obtain 



d_ 
di 



where 



Voru(t)J = VoVu(t) = Voru(t) • L^{t), 
Lu(i) = Vu(t)vu(t). 



(21) 



Taking into account that 



d r„ , ,1-1 

5^[v„r„(i)] 

we find from Eq. (21) that 

d_ 
di 



Voru(t) 



-1 d_ 
' di 



Voru(i) • Voru(t) 



Voru(t) =-K{t)- Vor,(t) 



(22) 



It follows from Eq. (13) that the derivative of the Cauchy deformation tensor C{t) reads 



dC 



where 



(i) = Vov(i)- Vor(i) +Vor(i)- Vov(i) 



(23) 



9 



is the velocity vector. Taking into account that 

Vov(i) = Vor(i) • V(i)v(i), 
where V(i) is the gradient operator in the deformed configuration, we arrive at the formula 

^(i) = 2Vor(t)-D(t)-[Vor(t)] , (24) 

where 



is the velocity gradient and 



L{t) = V(t)v(t) 



D(t) = i[L(t) + LT(t) 



is the rate-of-strain tensor. Equations (16), (20), (22) and (24) result in 



~~dt 



it) 



-K{t)-C,{t)-C,{t)-Ll{t) 



+2 



Voru(t) • Vor(t) • D(i) • Vor(i) • Wor^{t) 



(25) 



Bearing in mind that 



^^"(t) = -c-(t).^(t).c-(t). 



dt 



dt 



we find from Eqs. (16) and (25) that 



dt 

- 2[Voru(i)]^ • C-\t) ■ Vor{t) ■ D(i) • [Vor(i)]^ • C-\t) ■ Vor^{t). (26) 
The first principal invariant of the tensor Ce{t) is given by 

/l(Ce) = Ce : I, 

where the colon stands for convolution. Differentiating this equality with respect to time 
and using Eq. (25), we obtain 



dli 
'dt 



(Ce(t)) 



-L,(t):Ce(t)-Ce(t):L;[(t) 

+2[(Vor(i))^ • (Vor,(i))"^ • (Vor,(i))"' • Vor(i)] : D(i). 



Combining this equality with Eq. (14) and introducing the rate-of-strain tensor for sliding 
of junctions 

B,{t)^UK{t)+Ll{t)l 



we arrive at the formula 
dh 



dt 



(Ce(t)) = 2[(Vor(i)) • C-\t) ■ Vor(t)] : B{t) - 2Ce(t) : D,(t). 
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(27) 



It follows from Eqs. (19) and (27) that 

dh 



dt 



(Ce(t)) = 2[Be(t) : Bit) - Ce(i) : D,(t) 



(28) 



The Cauchy deformation tensors for transitions from the reference state to the deformed 
state and from the reference state to the unloaded state are assumed to obey the incom- 
pressibility condition, which implies that 

/3(Ce) = 1. 

It follows from this equality that 

/2(Ce(t)) = h{C;\t)) = C:\t) : I. (29) 
We differentiate Eq. (29) with respect to time, use Eqs. (13), (14) and (26), and obtain 
dh 



dt 



(Q 



C:\t):Uit) + Llit):C;\t)-2B{t): 

(Vor(i))^ • C-\t) ■ Voru(t) • (Voru(i))^ • C'^t) ■ Vor(t) 
2C;\t) : D,(t) - 2[(Vor(t))"' • C,(t) • (Vor(t 



: D(t). 



This equahty together with Eq. (19) results in the formula 
dh 



dt 



(Ce(i)) = -2[b;' : B{t) - C;\t) : D,(t) 



(30) 



Equations (28) and (30) determine the derivatives of the principal invariants of the tensor 

Ce{t). 

Our goal now is to calculate the derivative of an arbitrary smooth function, (j) = 
(f){h,h), of the principal invariants, h, of the Cauchy deformation tensor, Ce{t), with 
respect to time. According to the chain rule for differentiation. 



with 



dt ^'^ ' ^ dt ^' ^ ' ^ dt 



(31) 



Substitution of expressions (28) and (30) into Eq. (31) imphes that 

d(j) 



^(/l(Ce(t)),/2(Ce(t 



= 2[0,iBe(t)-0,2B;^(t)J :D(i) 
-2UiCe(i)-0,2C;^(t)l :D,(i). 



(32) 



Equation (32) generalizes the conventional formula for the derivative of a function of 
the principal invariants, hi of the Cauchy deformation tensor C(i), 



dt 



(/i(C(i)),72(C(i)))=2[0,iB(i)-0,2B-^(t) 
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D(t). 



(33) 



In what follows, we analyze deformation of meso- domains with various potential energies, 
cu, where junctions start to slip with respect to the bulk material at various times, r, 
while they moved together with the bulk medium before the instant r. This means that 
the tensors Be and Ce become functions of three variables, t, r and cu, that satisfy the 
initial conditions 



Be(i,T,a;)^_ =B(r), C,{t,T,uj) 



t=T 



t=T 



C(r). 



(34) 



It follows from Eq. (32) that the derivative of an arbitrary function of the first two 
principal invariants of the Cauchy deformation tensor Ce(i,T) reads 



dt 



(/i(Ce(t,r,a;)),/2(Ce(t,r,c^)) 
= 2 



(/.,iBe(t,T,u;) - (t>,2Q-\t,T,u;)\ : D(i) 

0,lCe(t, T, u) - (f),2C~^{t, T, u) : Du(i, T, u). 



(35) 



Equations (33) and (35) will be used in the next section to calculate the derivative of the 
mechanical energy of a filled elastomer with respect to time. 



4 Strain energy of an elastomer 

We confine ourselves to loading processes that increase the concentration of active mcso- 
domains with time. This implies that any active region can be entirely characterized by 
its potential energy, cu, and the instant, r, when it became active. We set r = for 
meso-domains active in the stress-free state. 

The conventional affinity hypothesis is adopted for deformation of passive meso- 
domains. This assumption disregards thermal oscillations of junctions and postulates 
that the deformation gradient for the motion of junctions at the micro- level coincides 
with that for the movement of appropriate points of an elastomer at the macro-level 
(Treloar, 1975). 

A passive region is treated as an isotropic incompressible medium whose mechanical 
energy per strand, w, depends on the current temperature, T (the strong effect of temper- 
ature reflects the entropic nature of the strain energy for rubbery polymers), the measure 
of damage, a, and the first two principal invariants, Ik, of the right Cauchy deformation 
tensor C, 

w = w{a{t),T{t), h{C{t)), HCit))). (36) 

At an arbitrary instant t > 0, the strain energy of a strand, wo, belonging to a meso- 
domain with potential energy u which became active at instant r G [0,t] is assumed to 
be a function of the current temperature, T, and the first two principal invariants of the 
right Cauchy deformation tensor Ce, 

^^;o = w^o(T(^),/l(Ce(t,r,o;)),/2(Ce(^,r,o;))). (37) 

Equations (36) and (37) are based on the conventional hypothesis that the excluded- 
volume effect and other multi-chain effects are screened for an individual chain by sur- 
rounding macromolecules and they may be accounted for in terms of the incompressibility 
condition (Everaers, 1998). 
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To simplify the analysis, the strain energy of a strand is assumed to be independent of 
the potential energy, a;, of a meso-region to which this strand belongs (the latter means 
that the functions w and wq are treated as average mechanical energies of strands). The 
mechanical energy of a strand in an active meso-region, wq, differs, however, from the 
strain energy, w, of a strand in a passive domain. The mechanical energy of a passive 
meso-domain with a potential energy uj and a damage measure a reads 

w{a, T, h, h) = w^{T, h, h) + aA«;(r, A, h), (38) 

where Aw is the increment of strain energy per strand in a passive meso-domain (with 
respect to that in an active region) caused by inter-chain interactions (van der Waals 
forces). Equation (38) implies that the strain energy density of an elastomer is continuous 
at the point of transition of a passive meso-domain into the active state (when the damage 
parameter, a, vanishes). 

Summing the mechanical energies of strands belonging to initial active meso-regions 
with various potential energies, a;, to passive meso-domains with various energies, a;, and 
damage measures, a, as well as to meso-domains that become active at various instants, 
T e [0, i], we find the strain energy density per unit mass of an elastomer 

W{t) = dujy^ Ep{t,a,uj)w{a,T{t),h{C{t)),l2{C{t)))da 
+X,{0, uj)wo{T{t), /i(Ce(t, 0, cu)), hiC.it, 0, cu))) 
+ 7(r, Lj)wo{T{t), /i(Ce(t, r, a;)), hiC.it, r, u;)))dT . 

Substitution of Eqs. (36) to (38) into this equality results in 

W{t) - iVp(t)ti;o(T(t),/i(C(i)),/2(C(i))) +Z(i)Aw(r(i),/i(C(i)),/2(C(i))) 
+ dw X^{0, u)wo{T{t), /i(Ce(t, 0, u)), hiCeit, 0, u))) 

+ 1^ 7(r, uj)wo{T{t), /i(Ce(t, r, ^)), hiC^it, r, w)))(ir] , (39) 
where the function Np{t) is given by Eqs. (6) and (7) and 

roc roo 

Z(t) — / ada / Ep(t,a,uj)dw. (40) 
Jo Jo 

Our purpose now is to calculate the derivative of the function W{t) with respect to time. 
It follows from Eqs. (7), (33) to (35) and (39) that 

^{t) = Jit)^it) + 2A(t) : B{t) - Y,{t) - Y,{t), (41) 

where 

At) = iVp(t)^(T(t),/i(C(t)),/2(C(t))) +Z(t)^(T(t),/i(C(t)),/2(C(i))) 
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roo 

+ duj 
Jo 



dWn 



X,{0, u;)-^(T{t),h{C,{t, 0, u)), /2(Ce(i, 0, u))) 



+ 1 7(r,a;)^(T(t),/i(Ce(t,r,u;)),/2(Ce(t,r,u;)))(ir], 
N^{t) {wo,i{T(t), h(C(t)), h(C{t)))B(t) 
-wo,2{m,h{C{t)),h{C{t)))B-\t)^ 
+Z(t) (a«;,i {T(t), h(C(t)), h(C(t)))B(t) 
-Aw,2(T{t),h{C{t)),h{C{t)))B-\t)^ 

+ duj Xa(0, uj) (wo,i {T{t), /i(Ce(i, 0, ou)) , hiC.it, 0, uj)))B,{t, 0, ou) 
-wo,2 {m, /i(Ce(t, 0, u;)), /2(Ce(t, 0, uj)))B;\t, 0, u;)) 
+ ^* 7(r, cu) (^u'o,i(T(t), /i(Ce(t, r, to)), hiC^it, r, Lo)))B,{t, r, u) 
-Wo,2 {T{t), /i(Ce(t, T, u)), /2(Ce(i, T, u;)))B,-i(i, T, u;)) drj , 



(it 
dZ 



^t)wo{T{t),h{C{t)),h{C{t))) 



^^{t)Aw{T{t),h{C{t)),h{C{t))) 

-fit, u;)wo{T{t), /i(Ce(t, t, a;)), /2(Ce(t, t, u;)))du;, 



Y2(t) = 2 / 

JO 



X,(0,u;)A„(i,0,u;) :D„(f,0,u;) 



+ / 7(r,a;)Au(t,r,a;) : Du(t, r, a;)c?r 



dcu. 



(42) 



The function A„ reads 



Au(t,r,a;) = u^o,i /i(Ce(t, r, cu)), /2(Ce(t, r, cc;)))Ce(t, r, cu) 

-w;o,2 (T(t) , /i (Ce(t , r, cu)) , /2 (Ce (t , r, uj)))C;\t, r, uj) . 
Equations (7), (9), (34) and (42) imply that 

Yiit) = -^{t)Aw(T{t),h{C{t)),UC{t))), 



(43) 



(44) 



which means that the function Yi is non-negative for an arbitrary program of loading with 
an increasing number of active meso-domains. 

5 Constitutive equations 

We confine ourselves to quasi-static loadings with finite strains, when mechanically in- 
duced changes in temperature, AT, are rather weak. This imphes that the effect of 
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temperature on the specific heat, c, may be disregarded. The following expression is 
adopted for the free energy (per unit mass) of a filled elastomer: 

■^^■^o + {c-So){T-To)-cT\n^ + W, (45) 

-to 

where and 5*0 are the free energy and the entropy per unit mass in the stress-free state 
at the reference temperature Tq. The second and third terms on the right-hand side of 
Eq. (45) characterize the free energy of thermal motion of chains. 

For an incompressible medium, the Clausius-Duhem inequality reads (Haupt, 2000) 

where p is mass density, S' is the deviatoric component of the Cauchy stress tensor S, q 
is the heat flux vector, S is the entropy, and Q is the internal dissipation per unit mass. 
Substitution of Eqs. (41) and (45) into Eq. (46) results in 



Q = _(^_^o_cln- + j)— + -(s'-2pA) :D 
1 



+n + i^2-^q-VT>0. (47) 

The function Yi describes the energy dissipation driven by transition of passive meso- 
domains into the active state, Y2 characterizes the entropy production induced by slippage 
of junctions with respect to the bulk material, whereas the last term in Eq. (47) is 
responsible for internal dissipation caused by thermal conductivity. 

Because Eq. (47) is to be satisfied for an arbitrary loading, the expressions in brackets 
vanish. This implies the formula for the entropy per unit mass 

S{t)^S^ + c\n^-J{t) 

and the constitutive equation 

S(i) = -P(i)I + 2pA(i), 

where P is pressure. Substitution of expression (42) into this equality implies the stress- 
strain relation 



S(i) = -P(i)I + 2p|7Vp(i)(^^o_,(r(i),/i(C(i)),/2(C(i)))B(i) 

+Z{t) (^Aw,,{T{t), h{C{t)), l2{C{t)))B{t) 
-A^,2(T(t),/i(C(t)),/2(C(t)))B-i(t)) 

X,{0, u) (wo,i {T{t), 7i(Ce(i, 0, u)), hiC.it, 0, u;)))B,{t, 0, u) 
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-Wo,2{T{t), /l(Ce(t, 0, LU)), /2(Ce(t, 0, Uj)))B;\t, 0, Uj) 

+ 7(r, u) (wo,i (T(t), /i(Ce(t, r, u)), h{C,{t, r, uj)))B,{t, r, cu) 
-«;o,2(r(i),/i(Ce(i,T,a;)),/2(Ce(i,T,u;)))B,-i(i,T,u;))ciT (48) 

We adopt the Fourier law for the heat flux vector q, 

q = -kVT, 

where k > is thermal diffusivity. It follows from this equation that the last term on the 
right-hand side of Eq. (47) in non-negative. 

The incompressibility condition implies that the first invariant of the rate-of-strain 
tensor Du vanishes, which means that 

A, : = : D„ 

where the prime stands for the deviatoric component of a tensor. By analogy with Eq. 
(1), we postulate that the rate-of-strain tensor for sliding of junctions, D^, is proportional 
to the tensor A(j, 

pA;(i,T,a;)=77D„(i,T,a;), (49) 

where > is a viscosity. Equation (49) provides the flow rule for slippage of junctions 
with respect to the bulk material, which ensures that the term A^ : Du is non-negative. 
It follows from this result and Eq. (42) that 1^2 > for an arbitrary program of strain- 
ing. Because the function Yi is non-negative as well, governing equations (48) and (49) 
guarantee that the dissipation inequality (47) is fulfilled. 
It follows from Eqs. (2), (4), (42) and (49) that 

pro 



Du(i,T, uj) 



G 



■ ex.p[—uj ) 



wo,i(T{t), h{C,{t, T, ou)), /2(Ce(t, T, u;)))c',{t, T, iu) 



-Wo,2(T(t), /i(Ce(t, r, a;)), hiC.it, r, uj))) {C;\t, r, u))' 



(50) 



Formulas (48) and (50) entirely determine the evolution of an ensemble of mcso-regions 
in a filled elastomer. These equations may be noticeably simplified provided that the 
conventional theory of rubber elasticity is applied to the description of the mechanical 
energy in active meso-domains. Setting 

pwo = G(/i-3), (51) 

which corresponds to the strain energy density of a neo-Hookean material, we find that 



E(i) = -P{t)I + 2G\N^it)B{t) + 



X,{Q,uj)B,{t,0,cv) 



+ [ J{T,UJ)Be{t,T,UJ)d^ 

Jo 



+2pZ{t) 



Aw,,{T{t),h{C{t)),h{C{t)))B{t) 



-Aw,2{T{t),h{C{t)),l2{C{t)))B-\t) 



(52) 
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and 



Du(t, T, uj) = To cxp(-cj)C^(t, T, u). 



(53) 



Constitutive equations (52) and (53) are applied in the next section to analyze stresses in 
a bar under tension. 



6 Uniaxial tension of a specimen 

Points of the bar refer to Cartesian coordinates {Xi} in the stress- free state and to Carte- 
sian coordinates {xi} in the deformed state, (i = 1,2,3). Uniaxial tension of an incom- 
pressible medium is described by the formulas 



Xi^k{t)Xi, X2^k ■^{t)X2, X3^k 2(^)^3, 



(54) 



where k = k{t) is the extension ratio. We assume that transition from the reference state 
to the unloaded state is determined by equations similar to Eq. (54), 



6 = ^u(^,T,a;)Xi, 5 = /cu ^(i,T,u;)X2, ^3 = /cu ' (i, r, u;)X3, 



(55) 



where {^j} are Cartesian coordinates in the unloaded configuration and k-a{t,T,u}) is a 
function to be found. The Cauchy deformation tensors B(i) and C(t) read 



B(i) = C{t) = k''{t)e^e^ + k-\t){e2e2 + 6363), 



(56) 



where are base vectors of the frame {Xi\. The Cauchy deformation tensors Be(t,r, a;) 
and Ce(t,T, a;) are given by 



Be(i, T, U) = Ce(i, T, u) = { } ) CiBi + 

V/Cu(t,r,u;)/ 

which implies that 



kn(t,T,Uj) 

k{t) 



(6262 + 6363), (57) 



C;(i,T,^) 



k{t) 



ku{t,T,u;) 



/Cu(t, r, uj) 

m 



6161 - -(6262 + 6363, 



(58) 



It follows from Eq. (55) that 
Du(t,T,u;) 



1 



T77 ~ ^(6262 + 6363) 

fcu(t, r, cu) ot L 2 



(59) 



Substitution of expressions (58) and (59) into the flow rule (53) implies that the function 
A;u(t, T, a;) satisfies the differential equation 



dku 



ku{t,T,uj) dt 
with the initial condition 



(t, r, a;) = r*exp(— a;) 





ku{t,T,uy 


\k^{t,T, uj), 


1 k{t) I 
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t=T 



(60) 
(61) 



The coefficient in Eq. (60) reads 



3 



It follows from Eqs. (52), (56) and (57) that the non-zero components, Ejt, of the Cauchy 
stress tensor 

S = SieiGi + ^2(6262 + 6363) 

are given by 

Ei(i) = -P(t) + 2pZ{t) [Aw iA;^(i) - Aw 2A;-^(i)] + 2Gk'^{t)^N^{t) 



-\t)[N,{t) 



E2(i) = -P{t) + 2pZ{t) Aw^ik-^t) - Aw^2k{t) 



+ 2Gk 



+ 



roo p rt -, ^ 

Xa(0,a;)A;u(i,0,a;) + -f{T,u;)k^{t,T,u;)dT duj, 



where the arguments of the function Aw are omitted. Excluding the unknown pressure 
P{t) from the expressions for Ej^ with the help of the boundary condition 



S2(t) = 0, 



we find the longitudinal stress 



Ei(t) 



GNp{t) + pZ{t)[Aw-^ + k-\t)Aw^2) {k^{t) - k-\t)) 



+2Gf [X.(0,.)((- 



k{t) 



(t,0,o^) 



k^{t,0,uj) 



kit) 



+ r7(r,.)((-^) 

Jo \^ku{t,T,UJ)^ 



dr 



dio. 



(62) 



Equations (60) to (62) determine the stress in a specimen for an arbitrary program of 
straining. To compare results of numerical simulation with experimental data, we focus 
on a tensile relaxation test with 



m = 



1, t<0, 
A, t>0, 



(63) 



where A > 1 is a constant. 

For the loading program (63), Aw 1 and Au'^2 become functions of A only, 

Awfe = Awfc(A) (A; = 1,2). 

The concentration of active meso-regions is assumed to be altered at the instant t — 
only, which implies that the function 7(^,0;) vanishes for any t > and the quantities A^p 
and Z are entirely determined by the elongation ratio A, 



Z{\). 
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The number, Xa(0, a;), of active meso-domains with a potential energy u) becomes a 
function of A, 

Xa(0,u;) = 7Va(A)p(A,u;), 

where Ng_{X) is the number of active regions at stretching to the elongation ratio A, and 
p(A, cu) is the distribution function for active meso-domains that satisfies the condition 

/ p(X,u})du} ^1. (64) 
Jo 

For the relaxation test (63), the quantity ku{t, 0,w) also becomes a function of A, 

k-a{t, 0, cu) — h{t, A, cu). 
Using Eqs. (10), (11) and introducing the notation, 

V'i(A) = 2 G'A^ + pZ(A)(Awi(A) + A-^Aw2(A)) (a^-A-^), 
^2(A) = 2GN,{\)[x'-\-'), 

we find from Eqs. (62) and (63) that 

Ei(i,A) = Vi(A)-V'2(A) / l-f{t,X,cu)p{X,cu)dcu, 

Jo L J 

where 



(65) 
(66) 

^h{t,x,uj)^ X y ^^^^ 

It follows from Eqs. (60), (61) and (63) that the function h{t,X,uj) obeys the differential 
equation 



dh 
dt 



(t. A, tu) = r^.exp(— cu) 





^2 h{t,X,uj)' 


\h{t, X,uj)/ 


1 X \ 



h{t,X,Lu), h{0,X,uj) = l. (68) 



Equations (64) and (66) imply that the dimensionless ratio 

Si(t,A) 



is given by 



where 



R(t, A) = 1 - A(X) / 1 - f{t, A, u) p(A, u)du, 

Jo L ^ 



To fit experimental data, we employ the quasi-Gaussian distribution function 

(a; - {u;{X))or 



p{X,uj) =po(A)exp 



2a^{X) 



(69) 
(70) 

(71) 
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where (a;)o and a are adjustable parameters and po is determined by condition (64). 
Substitution of expression (71) into Eq. (69) results in 



R{t,X)^l-A{X)po{X) / l-f{t,X,u;) 



(^-(^(A))o)^ 1 
2a2(A) . 

Introducing the new variable, z — cu — cu*, where a;* = In F*, we obtain 

{z-HxW^ 



duj. 



/oo _ 
l-f{t,X,z) 



dz, 



2(72 (A) 

where (a;) = (a;)o — UJ*. Assuming the relative width of the quasi-Gaussian distribution, 



c 



(72) 



to be small compared to {u) (this hypothesis will be discussed later), we neglect the 
integral over the interval [—a;*, 0] and find that 



R{t,X)^l-A{X)po{X) / l-f{t,X,z) 



2(t2(A) 



dz, 



(73) 



where the function f{t,X,z) is determined by Eq. (67). According to Eq. (68), the 
function h{t, A, z) satisfies the differential equation 



dh . ^ . . . 

— [t,X,z) = exp(-2;) 



\i ' ^ 


|2 h{t,X,z)^ 


\h{t, A, z)y 





h{t,X,z), h{0,X,z)^l. (74) 



Equations (73) and (74) are determined by 3 dimensionless parameters: 

1. the ratio. A, of the relaxing stress to the total stress at the beginning of the relaxation 
test; 

2. the average potential energy, (a;), for sliding of junctions in active in meso-domains; 

3. the standard deviation of the potential energy, a. 



7 Experimental procedure 

Three series of uniaxial relaxation tests were performed at room temperature. Dumb- 
bell specimens were provided by TARRC laboratories (Hertford, UK) and were used as 
received. The tests were executed on specimens of the compound EDS-16 on the base 
of natural rubber reinforced by high abrasion furnace black. The compound contains 45 
phr of carbon black with the mean diameter of filler particles of about 30 nm (Jager and 
McQueen, 2001). A detailed formulation of the compound is presented in Table 1. 

Experiments were carried out by using a testing machine designed at the Institute of 
Physics (Vienna, Austria) and equipped with a video-controlled system. To measure the 
longitudinal strain, two reflection lines were drawn in the central part of each specimen 
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before loading (with the distance 7 mm between them). Changes in the distance between 
these hnes were measured by a video-extensometer (with the accuracy of about 1 %). 
The tensile force was measured by using a standard loading cell. The nominal stress was 
determined as the ratio of the axial force to the cross-sectional area of a specimen in the 
stress-free state, 4 mm x 1 mm. 

To analyze the effects of pre-loading and thermal recovery, the following testing pro- 
cedure was chosen: 

1. Any specimen was loaded with the strain rate 6+ = 3.0 ■ 10^^ s^^ up to a given 
elongation ratio A, which was preserved constant during the relaxation test (tr = 1 
hour), and unloaded with the strain rate e_ = — e+. 

2. The specimen was annealed for tg, = 24 hours at room temperature. 

3. The specimen was loaded with the strain rate 6+ up to the elongation ratio Amax — 
4.0 and unloaded with the strain rate e_. 

4. The specimen was annealed for tg, at ambient temperature. 

5. The specimen was loaded with the strain rate e+ up to the same elongation ratio, 
A, as in the first test. This elongation ratio was preserved constant during the 
relaxation test, and, afterwards, the specimen was unloaded with the strain rate e_. 

6. The specimen was recovered in an oven for 24 hours at the constant temperature 
T — 100 °C and cooled down by air to room temperature for 4 hours. 

7. The specimen was loaded with the strain rate e+ up to the same elongation ratio as 
in the first test. This elongation ratio was preserved constant during the relaxation 
test, and, afterwards, the specimen was unloaded with the strain rate e_. 

The relaxation tests were carried out at four elongation ratios: Ai = 2.0, A2 = 2.5, 
A3 = 3.0, and A4 = 3.5. 

The longitudinal stress a was measured as a function of time t (the initial instant 
t = corresponds to the beginning of the relaxation test). For any elongation ratio, A, 
the dimensionless ratio R{t, A) is depicted versus the logarithm of time (log = log^g) 
Figures 1 to 3. 

These figures evidence that the longitudinal stress is noticeable decreased during the 
relaxation time, t,. — 1 h. For virgin specimens, this decrease equals 19.7 % at the 
elongation ratio Ai and it reaches 27.9 % at A4. 

After pre-loading to Amax and subsequent unloading, the time-dependent decrease in 
stress is weakened. The portion of the longitudinal stress relaxed during is estimated 
as 14.4 % at the elongation ratio Ai and as 20.0 % at A4. 

Thermal recovery of specimens at the elevated temperature T for = 24 h returns 
the specimens to their initial state. The decrease in the longitudinal stress during = 1 
h is estimated as 19.7 % at the smallest elongation ratio, Ai, and it grows up to 27.0 % 
at straining with A4. 

These changes in the time-dependent response of the CB filled natural rubber are 
associated with alternation of the internal structure of the compound at the micro-level: 
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1. the concentration of active domains per unit mass grows with an increase in the 
elongation ratio, 

2. it is noticeably reduced under pre-loading and subsequent unloading, 

3. the concentration of active meso- regions (practically) returns to its initial value after 
thermal recovery. 

To verify these assertions, we determine adjustable parameters in Eqs. (73) and (74) by 
matching experimental data and analyze the effect of the longitudinal ratio. A, on the 
coefficients A, {lu) and a. 



8 Comparison with experimental data 

Given an elongation ratio. A, the fitting procedure consists in the following. We fix 
intervals [0, (a;)inax] and [0, (Jmax] where the "best-fit" parameters (cu) and a are assumed 
to be found and divide these intervals into J subintervals by points {u})j — jA^^ and 
aj = jA„ with j = 1, . . . , J, = (c<j)inax/ J and A^- = cXmax/ J- For any pair of parameters 
{u!)i and (Tj, we find the constant po from condition (64), where the integral is evaluated 
numerically. Afterwards, we calculate the integral in Eq. (73) and determine the pre- 
factor A — A{i,j) [which ensures the best fit of the experimental curve Rexp{t, A)] by the 
least-squares method. As a measure of deviations between observations and results of 
numerical analysis, we chose the function 

^ihj) = Zl[^exp(^fe, A) - -Rnum(^fc, A)]^ (75) 
tk 

where the sum is calculated over all experimental points depicted in Figures 1 to 3 and 
the function i?num(^! A) is determined by Eq. (73). Finally, the "best-fit" parameters {cu) 
and a are determined as those that minimize functional (75) on the set 

{{uj)i,aj (i,j = 1,..., J)}. 

To ensure high accuracy in the approximation of observations, after determining the 
"optimal" parameters {Lj)i and aj, the above procedure is repeated for the new intervals 
and [(7j_i, Uj+i]. In the numerical analysis, we set (a;)i„ax = 10.0, iTmax = 
10.0 and J = 20. The integrals in Eqs. (64) and (73) are calculated by using Simpson's 
method with 100 points and the step Az = 0.2. For any z, the ordinary differential 
equation (74) is solved by the Runge-Kutta method with the time step At = 0.05 s. 

To assess the effect of straining of the average potential energy for sliding, (a;), and 
the standard deviation of energies, a, we plot these quantities versus the first invariant, 
Ji, of the Cauchy deformation tensor C. It follows from Eq. (56) that in tensile relaxation 
tests, 

/i = A^ + 2A-^ 
Experimental data are approximated by the linear functions 

(uj) ^uJo + uji{li - 3), a ^ao + (7i(/i - 3), (76) 
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where the coefficients ujk and ak are determined by the least-squares technique. These 
parameters are hsted in Table 2. Figures 4 and 5 show that phenomenological relations 
(76) ensure good accuracy of fitting. 

To provide some physical basis for Eq. (76), we refer to Eq. (51) that implies that for 
a neo-Hookean medium, the difference /i — 3 may be treated a measure of the mechanical 
energy of active meso-regions. Equations (76) show that the parameters (a;) and a of the 
distribution function (71) for potential energies of active meso-domains are proportional 
to their average strain energy. 

Equations (67) and (68) imply that /(O, A, a;) = 1 and /(oo,A,a;) = 0. These equali- 
ties, together with Eqs. (64) and (66), yield 

Ei(0, A) = Vi(A), Ei(oo, A) = Vi(A) - M>^)- (77) 

We associate Ei(oo,A) with the equilibrium stress, Eeq(A), and AE(A) = Ei(0, A) — 
Ei(oo, A) with the relaxing stress. It follows from Eqs. (70) and (77) that the ratio 

reads 

a(A) — 



1 - A{\) 

Given an elongation ratio A, we find the dimensionless ratio a(A) from this equality, where 
the value of A{\) is determined by fitting observations, and plot a versus the first invariant, 
/i, of the Cauchy deformation tensor. Figure 6 demonstrates that the experimental data 
are correctly approximated by the linear relation 

a = oo + ai(/i — 3), (79) 

where the coefficients (found by the least-squares algorithm) are collected in Table 2. 

Figure 4 demonstrates that the average potential energy for sliding of junctions in- 
creases with the elongation ratio A. This may be explained by the fact that only meso- 
domains with low potential energies are active in the stress-free state, whereas regions with 
higher potential energies become involved in the sliding process provided that sufficiently 
large tension is applied to the specimen. The rate of growth of the average potential 
energy for sliding, (uj), with Ji is relatively small for virgin specimens and it noticeably 
grows after pre-loading. After thermal recovery, changes in the function (a;)(/i) become 
similar to those for a virgin sample. The loading-unloading procedure with Amax = 4.0 
results in a decrease of the average potential energy of the ensemble of active meso-regions 
for a stress-free specimen (the parameter ujq is reduced by 8 %). This observation may be 
explained by (at least) two reasons: 

1. damage of the rubbery compound during pre-loading "turns off' some active regions 
in a virgin sample that become passive under unloading; 

2. stretching of a specimen to a high elongation ratio and subsequent unloading increase 
the rate of slippage of junctions with respect to the bulk medium. According to Eq. 
(68), this mechanically induced acceleration of the sliding process is tantamount to 
a decrease in the potential energy of active meso-regions. 
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Figure 5 shows that the standard deviation of potential energies, cr, increases with the 
elongation ratio, A. This growth is rather modest for virgin and recovered samples, but 
it becomes stronger after pre-loading: the coefficient ui for the pre-loaded specimens 
exceeds that for the virgin specimens by a factor of 6. After stretching of specimens 
to Ainax and subsequent unloading, the quasi- Gaussian distribution of potential energies 
becomes substantially sharper: the standard deviation of energies for a stress-free rubbery 
compound, ao, is decreased by 40 % compared to the virgin material. This may be 
explained by the fact that unloading induces transformation of active meso-domains with 
relatively high potential energies (which characterize a "tail" of the Gaussian function with 
large values of uj) into the passive state. Thermal recovery at an elevated temperature 
"activates" a part of these regions (the corresponding curve in Figure 5 is located higher 
than the curve for samples suffered pre-loading), but no total recovery is observed: the 
graph of the function a{Ii) for the recovered compound does not reach that for the virgin 
medium. 

To evaluate the dimensionless width of the quasi-Gaussian distribution function (71), 
we calculate the ratio ( with the help of Eq. (72) and plot it versus the first invariant, Ji, 
of the Cauchy deformation tensor C in Figure 7. Experimental data are approximated 
by the linear dependence 

C = Co + (1^1 -3), (80) 

where the constants Co and Ci are determined by the least-squares algorithm. Figure 7 
demonstrates fair agreement between the observations and the predictions of the phe- 
nomenological equation (80) with the parameters Co a-^d Ci listed in Table 2. 

According to Figure 7, the dimensionless ratio ( is located within the interval between 
0.30 and 0.45 for all specimens, which imphes that the width of the distribution function 
(71) is relatively small. The latter may be thought of as a confirmation of the assumptions 
employed in the derivation of Eq. (73). The quantity ( slightly decreases with the 
elongation ratio A for virgin and recovered specimens, and it weakly increases with A for 
pre-loaded samples. The maximal width of the quasi-Gaussian distribution is observed 
for the virgin material. The distribution of potential energies for sliding of junctions in 
active meso- regions, p(a;), becomes sharper after pre-loading, whereas thermal recovery 
results in widening of this distribution. 

The ratio of the relaxing stress to the equilibrium stress, a, is plotted in Figure 6 as a 
function of the extension ratio, A. The parameter a is strongly influenced by stretching: 
for example, it grows by 45 % for virgin specimens when A increases from 2.0 to 3.5. 
Equations (65), (77) and (78) imply that 

,(^^ ^ N 

N^{X) + pG-^Zi\)[Aw,,{X) + X-^Aw4X)y ^ ' 

Because the term in brackets in Eq. (81) [the extra stress caused by inter- molecular 
interactions in passive meso-domains] is an increasing function of A, the growth of a(A) is 
associated with the following phenomena occurring simultaneously: 

1. micro-damage of passive meso-regions which results in a decrease of Z{X), 
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2. transformation of passive domains into the active state. The latter is equivalent to 
an increase in iVa(A) and a corresponding decrease in Np{\) that is connected with 
Na{X) by the conservation law (10). 

To assess the effect of these two sources for changes in a(A), constitutive relations are 
necessary for damage of inter-chain links in passive meso-regions. A model for rupture of 
van der Waals links under loading will be the subject of a subsequent publication. 

Stretching of samples to Amax and subsequent unloading result in a pronounced de- 
crease in a that is partially recovered after annealing at an elevated temperature. The 
severe decrease in a (the ratio for the relaxing stress to the equilibrium stress in pre- 
loaded specimens at relatively small strains, oq, equals about 50 % of its value for virgin 
samples) reflects the decay of the relaxation process evidenced by comparison of curves 
plotted in Figures 1 and 2. This decrease in a is assumed to be closely connected with the 
Mullins effect observed in loading-unloading tests on particle-reinforced elastomers (see, 
e.g., Bouche, 1960, 1961; Govindjee and Simo, 1992; Bergstrom and Boyce, 1998; Miehe 
and Keck, 2000; Wu and Liechti, 2000, to mention a few). 

It is worth noting a particular case of the model where the function Sp is independent 
of time (or the current strain for a time-dependent program of loading). This implies that 
the parameters A^aj Ap and Z remain constants, and 7 vanishes. Under these assumptions, 
the stress-strain relations (52) and (53) are reduced to the constitutive equations for a 
rheological model consisting of a nonlinear spring (which is characterized by the strain 
energy Aw of inter-chain interactions) connected in parallel with an infinite array of the 
Maxwell elements [nco-Hookcan springs linked in series with nonlinear dashpots whose 
response is described by Eq. (53)]. A variant of the latter model with a finite number 
of Maxwell's elements was widely used in the analysis of the viscoelastic and viscoplastic 
behavior of filler rubbers, see, e.g., Govindjee and Simo (1992), Kaliske and Rothert 
(1998), Miehe and Keck (2000), Haupt and Sedlan (2001). Fitting observations in the 
relaxation tests reveals that this simplified model does not capture noticeable changes 
in the parameters (uj), a and a under stretching. Our approach may be treated as an 
extension of the conventional rheological model, where the number of Maxwell's elements 
and their viscosities are affected by mechanical factors. 

9 Concluding remarks 

Constitutive equations have been derived for the time-dependent behavior of particle- 
reinforced elastomers at finite strains. A filled rubber is modelled as an ensemble of 
meso-domains where junctions between chains (chemical and physical cross-links and en- 
tanglements) are not fixed in the bulk medium, but they can slip with respect to it. A 
mean-field approximation is employed to develop stress-strain relations: a complicated 
micro-structure of an elastomer reinforced with aggregates of filler is replaced by equiva- 
lent networks of macromolecules in meso-regions with various potential energies for sliding 
of junctions. The distribution of potential energies of these meso-domains refiects 

1. impurities and local non- homogeneities in the spatial distribution of a cross-linker 
(which is typical of both unfilled and filled elastomers). 
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2. density fluctuations driven by segregation of sliort cliains to tlie interfaces between 
filler particles and the host medium, 

3. severe micro-stresses in the neighborhoods of filler clusters in a strained elastomer 
caused by the strong difference in the elastic moduli of filler and rubber. 

It is assumed that junctions slip with respect to the bulk material in active meso-domains, 
and the rate of sliding is determined by Boltzmann's formula. Sliding of junctions in 
passive meso-regions is prevented by surrounding macromolecules (by means of inter- 
chain links that are partially ruptured under loading). The concentration of links that 
forbid sliding is described by a measure of micro-damage, a. Activation of a passive 
meso-domain occurs when the van der Waals forces between strands vanish that prevent 
sliding of junctions. 

The mechanical energy of a filled elastomer is determined as the sum of strain ener- 
gies of chains in active and passive meso-domains and the energy of interaction between 
macromolecules in passive regions. Constitutive equations for a filled elastomer are de- 
rived by using the laws of thermodynamics (the study is confined to the loading processes 
in which the concentration of active meso-rcgions docs not decrease with time). These 
equations are applied for the analysis of stresses in an incompressible bar under uniaxial 
tension. 

To vahdate the model, three series of tensile relaxation tests were performed on carbon 
black filled natural rubber at elongations up to 350 %. In the first series, virgin specimens 
are used, in the second series, relaxation experiments were carried out on the samples 
after mechanical pre-loading, and in the third series, the time-dependent response was 
measured on the same specimens after recovery at an elevated temperature. Figures 1 to 
3 demonstrate fair agreement between the experimental data and the results of numerical 
simulation. 

To describe quantitatively changes in the internal structure of a particle-reinforced 
elastomer induced by pre-loading and thermal recovery, the distribution of potential en- 
ergies for sliding of junctions is approximated by a quasi-Gaussian function with two 
adjustable parameters. It is revealed that the average potential energy and the stan- 
dard deviation of energies grow with the elongation ratio A. After the loading-unloading 
procedure, these parameters are noticeably reduced. The latter is associated with rear- 
rangement of active meso-domains driven by rupture of filler clusters. Thermal recovery 
implies that the parameters of the quasi-Gaussian distribution attempt to return to their 
values for virgin specimens, but they do not reach their initial values after annealing for 
= 24 h at T = 100 °C. 
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Table 1: Chemical composition of compound EDS-16 (phr by weight) 



Formulation 


Content 


Natural rubber 


100 


Zinc oxide 


5 


Stearic acid 


2 


Carbon black (N330) 


45 


Process oil 


4.5 


Antidegradant(HPPD) 


3 


Antiozonant wax 


2 


Accelerator (CBS) 


0.6 


Sulphur 


2.5 



Table 2: Adjustable parameters of the model 



Specimens 








0-1 






Co 


Ci 


Virgin 


7.508 


0.109 


3.252 


0.017 


0.296 


0.028 


0.423 


-0.002 


Preloaded 


6.922 


0.182 


2.031 


0.103 


0.147 


0.018 


0.293 


0.006 


Recovered 


7.397 


0.134 


2.871 


0.019 


0.280 


0.025 


0.375 


-0.002 
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Figure 1: The dimensionless ratio R versus time t s in a relaxation test with the elongation 
ratio A (virgin specimens). Circles: experimental data. Solid lines: results of numerical 
simulation 
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Figure 2: The dimensionless ratio R versus time t s in a relaxation test with the elongation 
ratio A (the specimens after pre-loading). Circles: experimental data. Solid lines: results 
of numerical simulation 
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Figure 3: The dimensionless ratio R versus time t s in a relaxation test with the elongation 
ratio A (the specimens after recovery). Circles: experimental data. Solid lines: results of 
numerical simulation 
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Figure 4: The average potential energy for sliding of junctions (cu) versus the first invariant 
h of the Cauchy deformation tensor. Symbols: treatment of observations. Unfilled circles: 
virgin specimens; filled circles: the specimens after pre-loading; diamonds: the specimens 
after recovery. Solid lines: approximation of the experimental data by Eq. (76) 
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Figure 5: The standard deviation of potential energies a versus the first invariant l\ of the 
Cauchy deformation tensor. Symbols: treatment of observations. Unfilled circles: virgin 
specimens; filled circles: the specimens after pre-loading; diamonds: the specimens after 
recovery. Solid lines: approximation of the experimental data by Eq. (76) 
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Figure 6: The dimensionless parameter a versus the first invariant /i of the Cauchy de- 
formation tensor. Symbols: treatment of observations. Unfilled circles: virgin specimens; 
filled circles: the specimens after pre-loading; diamonds: the specimens after recovery. 
Solid lines: approximation of the experimental data by Eq. (79) 



36 



1.0 



0.0 




0.0 



/i-3 



10.0 



Figure 7: The dimensionless parameter C versus the first invariant Ii of the Cauchy de- 
formation tensor. Symbols: treatment of observations. Unfilled circles: virgin specimens; 
filled circles: the specimens after pre-loading; diamonds: the specimens after recovery. 
Solid lines: approximation of the experimental data by Eq. (80) 
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